Melting a copper cluster: Critical droplet theory 
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We simulate the melting of a 71 A diameter cluster of Cu. At low temperatures 
the crystal exhibits facets. With increasing temperatures the open facets pre-melt, 
the melted regions coalesce into a liquid envelope containing a crystalline nucleus, 
and the nucleus finally goes unstable to the supercooled liquid. Using critical droplet 
theory and experimental data for Cu, we explain the thermodynamics of the coexis- 
tence region. The width of the transition scales as (Number of particles) -1 / 4 . 
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We simulate the melting of an approximately spherical, 71 A diameter, 16727 atom 
cluster of Cu, using molecular dynamics with effective medium theory interactions. [|1],|2| By 
using a realistic interaction potential we are not only able to make specific experimental 
predictions, but we are also able to compare our results to analytical calculations using 
experimentally determined latent heats, specific heats, surface tensions, etc. We analyze 
our results in three contexts. First, we observe many of the effects predicted in the theory 
of equilibrium crystal shapes, || including liquid regions on open facets which can either 
be interpreted as pre-melting J| or "liquid lenses". |||6[] Second, we use critical droplet 
theory [F| to quantitatively explain the thermodynamics of the transition. By working at 
fixed energy (rather than fixed temperature) the first-order transition unfolds, revealing an 
interesting partially melted region in which adding more energy to the cluster reduces its 
mean temperature (the system has a negative specific heat). Third, we make contact to 
the formal theory of finite-size effects in first-order phase transitions. We show that the 
broadening we observe in the transition scales as the number of atoms A -1 / 4 , which for large 
iV dominates the other studied broadening mechanisms. || 

The molecular dynamics simulations of the properties of this Cu cluster are carried out 
using the Verlet algorithm for time-integration of the equation of motion. The simulations 
employ Andersen thermalization || for obtaining a given temperature or total energy of the 
system. The calculations were done on a massively parallel Connection Machine CM-200 
computer, using an algorithm where the face-centered cubic computational box was subdi- 
vided into a 16 3 grid. A computational processor is assigned to each grid point computing 
for the atoms in its associated grid sub-box. The interatomic interactions of the effective- 
medium theory (which extend roughly to fourth neighbors) are carried out by regular com- 
munication of atomic properties inside 5 3 sub-boxes surrounding each atom, assuring that 
all non-zero interactions are accounted for correctly, as in the serial-computer algorithm in 
ref. 0. These algorithms will be described in a succeeding paper. [[IT]] 

In the present study of the solid-liquid phase transition, it is natural to control the 
total energy parameter and measure the temperature as the average kinetic energy of the 
atoms. The fluctuation ((AT) 2 ) of the equilibrium temperature T at constant energy can 
be shown to be k B T 2 [2/ (3Nk B ) — 1/C], where A is the number of atoms, and C denotes 
the specific heat of the cluster as a whole. The measured fluctuations agree well with this 
formula except near the instability point where the critical slowing down (see below) makes 
sufficient statistics impractical to collect. 

The Cu cluster was initially a spherically truncated fee crystal, heated to 1282 K in 
order to melt about 3 layers of surface atoms. This cluster was both cooled and heated, 
respectively, at a sufficiently slow rate to maintain thermal equilibrium. Figure 1 shows 
the facet distribution developing after the cooling. The relative areas of these facets in 
equilibrium depend on their relative surface energies through the Wulff construction, || and 
are therefore quite sensitive to the model of interatomic interactions. The effective-medium 



theory has been shown to describe the various Cu surface energies with good accuracy. [11 



Figure 1. The faceted copper cluster at low temperature (522 K), with prominent (111) 
and (100) facets, clear areas of (110) facets, and terraces on the (111) facets. 

Figure 2 shows the gradual melting of the cluster as the energy is increased. Figure 2a 
shows pre-melted (110) and (100) facets, with the (110) regions appearing to form liquid 
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lenses similar to those predicted by Lowen. |J At slightly higher energies, the melted regions 
coalesce into a liquid envelope containing a crystalline nucleus. Figure 2b shows the crys- 
talline nucleus near the energy at which it becomes unstable. We also did a simulation series 
starting from a fully melted liquid cluster with a minuscule central fee seed. This formed a 
highly defective cluster with many stacking faults which slightly modified the pre-melting 
behavior, but left the remainder of the thermodynamics virtually unchanged. 

Figure 2. Two cross sections of the cluster, at energies (a) -3.065 and (b) -2.990 
eV/atom. The cutting planes shown are (a) (100) and (b) (111). The atomic positions 
are averages over the equilibrated part of the simulation. The liquid atoms are shown in 
black, the crystalline atoms in white. A histogram of position fluctuations was used to 
distinguish liquid from crystalline atoms. 

The qualitative behavior of the dynamics in this simulation is also of interest. In the case 
of the highly defective cluster, we observed one of the stacking faults relieve itself during 
the simulation, when a partial dislocation nucleated and passed through the crystal. The 
dynamics of defect creation and annihilation will bear further study. For both simulations 
series, the surface diffusion is rather rapid on the melted surfaces, but moderately slow 
on the crystalline surfaces on the time scales investigated. Finally, we observe significant 
critical slowing down of the equilibration between crystal and liquid near the instability 
point. Deep in the liquid and crystalline regions of temperature, the equilibration time- 
scales were essentially zero; at the energy shown in figure 2b, they slowed down to several 
tens of picoseconds (requiring several days of Connection Machine time for full equilibration). 

In figure 3a, we plot E(T), the total energy of the cluster versus its temperature T. For 
an infinite system at constant temperature, this plot would consist of two straight lines and 
an abrupt vertical jump at T c (ignoring the exponentially small, experimentally unobservable 
effects of droplet fluctuations). For the present cluster, this simple jump unfolds. The upper 
line represents the liquid state. Below around 1315 K the liquid becomes metastable, but 
the nucleation barrier to form a crystalline nucleus remains insurmountable until much lower 
temperatures. 

The straight portion of the lower curve at low temperatures represents the faceted crystal. 
The slope of the curve starts to increase well below T c , representing the latent heat of 
pre-melting for the various facets, starting with the open (110) surface. There have been 
experimental observations of superheating in metal clusters, [|TI| which could be due to 
pinning of the liquid lens boundaries to the facet edges. Superheating has also been 
observed in simulations [O] for smaller magic-number clusters, but it is likely that this 



superheating reflects the lack of steps or adsorbed atoms on their facets. We have observed 
no signs of superheating or latent heat jumps here. It is possible that the pinning barriers are 
small in our cluster, but may become observable for the much larger experimental clusters. 

We find that the crystalline nucleus detaches completely from the surface of the cluster 
below the temperature maximum in the lower curve. At higher energies, the temperature de- 
creases as the energy is increased (as observed in simulations of smaller clusters by Labastie 
and Whetten [14]]). This is a simple consequence of critical droplet theory. The free energy 
per atom of a cluster with a solid nucleus containing the fraction rj of the atoms in the cluster 
is approximately f(r/) = —L[{T c — T)/T c ]r] + 'yr] 2 ^ 3 . The first term represents the free energy 
difference between the supercooled liquid and the crystal and L is the latent heat per atom. 
The second term represents the solid-liquid interface tension and the coefficient 7 is propor- 
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tional to a suitable average of the interfacial free energy density over the equilibrium crystal 
shape. The forces from these two terms balance when (T c — T)/T c = (2/3) / L)^ 1 ^ 3 so the 
smaller the crystalline nucleus, the larger the undercooling needed to stabilize it. When en- 
ergy is added to the system, part of the crystal melts, absorbing the energy. The remainder 
of the crystal has a smaller radius of curvature, and the interface tension demands a larger 
undercooling: even more of the crystal melts leading to a net decrease in temperature. 

Figure 3. In (a) is shown the E(T) energy vs. temperature plot for the copper cluster 
simulation at constant energy. The black points represent the energies depicted in Fig. 2. 
In (b) is shown E(T) as derived from our simple model, Eq. p.2| . Here, the lower solid curve 
is the (partially) crystalline state, the upper solid curve is the liquid state, the dashed curve 
represents the unstable critical nucleus. The dotted horizontal line denotes the energy at 
which the entropies of the liquid and the mixed phases are equal. 

A rough estimate of the temperature range over which the transition is broadened can 
be found by estimating the undercooling at which the crystalline nucleus becomes unstable. 
Equating the temperature drop from melting a small region to the shift in the undercooling 
needed to stabilize the now smaller crystalline region yields an instability undercooling AT ~ 
(7 3 /L 2 ciVT c ) 1//4 ~ iV~ 1//4 , where here c represents a weighted average of the specific heats 
of the bulk liquid and solid. This broadened transition is only visible in constant-energy 
simulations; as iV — > oo, it is larger than the finite-size effects at constant temperature 
considered in previous work. || Using the more complete analysis below, we have tested the 
size dependence of both the crossover temperature (where the liquid entropy equals that of 
the cluster) and the instability temperature (where the cluster becomes absolutely unstable): 
both scale as iV _1//4 . 

We can turn this simple explanation into a quantitative calculation. First, because we 
work at constant energy, thermodynamics tells us that we must maximize the entropy per 
atom s (rather than minimize the free energy /): 

s(r), e s , e{) = i]s s (e s ) + (1 - r/)s z (e z ). (0.1) 

We here regard the entropy as a function of the solid fraction 77, the energy per atom in the 
solid nucleus e s , and the energy per atom in the liquid region e z . The entropy in the solid 
phase s s is given as function of energy by s s (e s ) = s c s + c s log[(e s — e c s )/(T c c s ) + 1] where 
we make the approximation that the specific heat c s is independent of temperature. The 
superscript c indicates values at the transition point. The entropy in the liquid phase can 
be similarly defined, using that the energy difference per atom between the two phases at 
the critical temperature is given by the latent heat L = el — e c s , and the entropy difference 
iB8f-sZ = L/T e . 

The energy per atom e we write as 

e(ij, e„ e { ) = V e s + (1 - vh + j^[R 2 Jiv + R 2 s (lsi + A 7 exp(-2(i? - R,) /£))], (0.2) 

with A7 = 7s,; — 7 s i — 7 to where j sv , 7^, and ji v denote the free energy densities of the 
solid-vapor, solid-liquid, and liquid-vapor interfaces, respectively. The radius of the cluster 
R and the radius of the solid nucleus R s can be expressed in terms of the densities and the 
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solid fraction 77. The last term in Eq. p.2| represents the interaction between the solid-liquid 
and the liquid-vapor interfaces and gives rise to premelting of a surface if A7 is positive. As 
usual, [jlj] the ambiguity of where to place the liquid-solid interfacial position is reflected in 
the breakup of the surface free energy into entropy and energy: we use the convention that 
attributes the free energy cost entirely to energy 

The experimental values of most of the parameters entering Eqs. p.l| and p.2| are known. 
n|[rT|] There is substantial uncertainty in only two: The solid-vapor and solid-liquid in- 



terfacial energies. The value for the solid-liquid interfacial energy 7^ we take as |17j 263 
mJ/m 2 . However, values as low as 177 mJ/m 2 have been found. The value for the 



solid-vapor interfacial energy ^ sv is poorly known. Here, we take the value from ref. ||17[| . 
A change in this parameter of only 1 percent will change the A7 by a factor of 2. This pa- 
rameter controls the onset of pre-melting, and hence the shape of E(T). With the present 
model (ignoring anisotropy as well as faceting and disorder on edges and vertices) and the 
given experimental values, we do not expect to reproduce in detail the pre-melting observed 
in the simulations. The exponential decay of the interaction between the solid-liquid and 
liquid-vapor interfaces follows from a Ginzburg-Landau analysis ]I7| of pre-melting of flat 
surfaces. The correlation length £ can be estimated from the solid-liquid interfacial free 
energy density using the Hansen- Verlet melting rule. [17 Maximizing the entropy, Eq. p.l| , 



with respect to e s and t\ at fixed energy e, Eq. pT2| , gives rise to the natural condition that 
the solid and liquid parts of the cluster must have the same temperature, and a further max- 
imization with respect to the solid fraction 77 leads to the E(T) curves shown in figure 3b. 
At constant temperature, the transitions are vertical on this plot: there are two metastable 
states separated by the critical nucleus. The pre-melted crystal becomes unstable when the 
E(T) curve has a vertical tangent. At constant energy, on the other hand, the crystalline 
nucleus surrounded by liquid becomes unstable when the tangent is horizontal. 

Comparing figures 3a and 3b, we find that the molecular dynamics simulations using 
effective-medium theory potentials are in surprisingly good agreement with our simple model 
that uses only experimental data for the specific heats, the latent heat, and the effects of 
surface tension. However, there is a fairly substantial temperature shift between the figures. 
The experimental transition temperature for bulk copper is 1356.2 K: the interfacial tensions 
depress the transition in the model cluster to about 1190 K as seen in figure 3b. On the 
other hand, the transition is seen at 1335 K in the simulations. This could in part be due to 
the large uncertainty in the experimental value for the crystal-liquid surface tension used in 
our simple model, but it is also conceivable that the effective-medium theory JT|fJ] transition 
temperature could be off by the necessary 145 K. Simulations for a range of cluster sizes 
could pinpoint these properties. 

We have shown that simple critical droplet theory, developed to study nucleation rates, 
provides a complete explanation for our copper cluster melting problem as studied through 
molecular dynamics simulations. We conclude by using the theory to calculate the nucleation 
rate for constant total energy. The dotted horizontal line in figure 3b shows the energy 
e C ross a t which the entropy of the liquid state equals that of the mixed state, which at 
this point has about 7200 atoms in the crystalline nucleus. Below e cross the liquid droplet 
is supercooled and metastable, as is the mixed phase above this energy. Near e cross , the 
system is in principle in a mixed state, where the crystalline nucleus appears and disappears 
with thermal fluctuations. The time needed to see a fluctuation, though, is huge. The 
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intersection of the dotted line with the middle, dashed curve represents a transition state 
with a smaller, crystalline nucleus of the same energy and consisting of approximately 1100 
atoms. However, this critical nucleus has an entropy lower by about 170k B , which is thus 
the entropy barrier between the two metastable states. The nucleation rate of the crystal 
from the supercooled liquid is consequently some microscopic prefactor times e~ 170 . We 
conclude that the nucleation rate is negligible both in simulations and experimentally until 
temperatures much lower than T c . 
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